##### Step 1: load data and packages #####



library(tidyverse)
library(lubridate)
library(data.table)
library(haven)
library(fixest)
library(magrittr)
library(ggplot2)
library(MASS)
library(plm)
library(lmtest)
library(broom)
library(xtable)

base_dir <- "..." # Manually change this home directory
data_dir <- paste0(base_dir, "data/")

options(scipen=999)

## function to more easily specify controls in regression formulas
fixest_form <- function(
  y,
  x, 
  fe="0",
  x_re = FALSE,
  df = NULL
  ) {

  if (x_re & !is.null(df)) {
    x <- x %>% map(~ grep(., colnames(df), value=T)) %>% reduce(c)
  }
  
  xpart <- paste(x, collapse = " + ")
  fepart <- paste(fe, collapse = " + ")



  y %>% 
    paste(paste(xpart, fepart, sep="|"), sep=" ~ ") %>%
    as.formula

}

chpos <- readRDS(paste0(data_dir, "cable_zip.rds"))
load(paste0(data_dir, "zip_demographics_2010.RData"))

census2010 %>% setDT
chpos <- census2010[chpos, on = .(zipcode)]

chpos[,zip:=as.integer(zipcode)]

zip_fips <- fread(paste0(data_dir, "DMA-zip.csv")) %>% 
  .[,.(fips=FIPS, zipcode = str_pad(ZIPCODE, pad="0", width=5), dmacode=`DMA CODE`)] %>% 
  unique(on="zipcode")
zip_fips %<>% mutate(
  fips = str_pad(fips, pad="0", width=5)
)
zip_fips %<>% mutate(
  state = as.factor(str_sub(fips, 1, 2))
)
n_distinct(zip_fips$state)
chpos %<>% left_join(zip_fips, by="zipcode")
chpos <- chpos[position_year==2009]

tpcontrib10 <-
  fread(paste0(data_dir, "tea-party-summ-all-time-varying-house.csv"))
nontpcontrib10 <-
  fread(paste0(data_dir, "non-tea-party-summ-all-time-varying-house.csv"))
contrib08 <-
  fread(paste0(data_dir, "summ-r-house-2008.csv"))
demcontrib <-
  readRDS(paste0(data_dir, "donations-congressional-summ.rds"))

setnames(tpcontrib10, old = "contributor_zipcode_cleaned_5", new="zip")
setnames(nontpcontrib10, old = "contributor_zipcode_cleaned_5", new="zip")
setnames(contrib08, old = "contributor_zipcode_cleaned_5", new="zip")
setnames(demcontrib, old = "contributor.zipcode.cleaned.5", new="zip")





##### Step 2: combine Nielson, DIME, and Census data #####



chpos <- tpcontrib10[chpos, on = .(zip)]
chpos <- nontpcontrib10[chpos, on = .(zip)]
chpos <- contrib08[chpos, on = .(zip)]

demcontrib %<>% dplyr::select(zip, amount_dem, n_dem)
demcontrib %>% setDT()
demcontrib$zip %<>% as.integer()
chpos <- demcontrib[chpos, on = .(zip)]

chpos %<>% as_tibble()

  chpos %<>%
    rename(amount_tea_party = amount_tea_party_house,
           amount_non_tea_party = amount_non_tea_party_house,
           n_donors_tea_party = n_donors_tea_party_house,
           n_donors_non_tea_party = n_donors_non_tea_party_house)

chpos %<>% setDT()

chpos[is.na(amount_tea_party) | amount_tea_party<0, amount_tea_party:=0]
chpos[is.na(amount_non_tea_party) | amount_non_tea_party<0, amount_non_tea_party:=0]
chpos[,tpfrac:= amount_tea_party / (amount_tea_party + amount_non_tea_party)]
chpos[is.nan(tpfrac), tpfrac:=NA]
chpos[is.na(n_donors_tea_party) | n_donors_tea_party<0, n_donors_tea_party:=0]
chpos[is.na(n_donors_non_tea_party) | n_donors_non_tea_party<0, amount_non_tea_party:=0]

chpos[is.na(amount_dem) | amount_dem<0, amount_dem:=0]
chpos[is.na(n_dem) | n_dem<0, n_dem:=0]

chpos %<>% mutate(
  log_amount_tea_party = log(amount_tea_party+1),
  log_amount_non_tea_party = log(amount_non_tea_party+1),
  log_n_donors_tea_party = log(n_donors_tea_party+1),
  log_n_donors_non_tea_party = log(n_donors_non_tea_party+1)
)

chpos %<>% mutate(
  state = case_when(
    zip == "29707" ~ "45",  # South Carolina
    zip == "33472" ~ "12",  # Florida
    zip == "33578" ~ "12",  # Florida
    zip == "33579" ~ "12",  # Florida
    zip == "60642" ~ "17",  # Illinois
    zip == "73012" ~ "40",  # Oklahoma
    zip == "73025" ~ "40",  # Oklahoma
    zip == "78633" ~ "48",  # Texas
    zip == "78665" ~ "48",  # Texas
    zip == "80023" ~ "08",  # Colorado
    zip == "84081" ~ "49",  # Utah
    zip == "85138" ~ "04",  # Arizona
    zip == "85142" ~ "04",  # Arizona
    zip == "85143" ~ "04",  # Arizona
    zip == "85286" ~ "04",  # Arizona
    zip == "85295" ~ "04",  # Arizona
    zip == "85298" ~ "04",  # Arizona
    zip == "85392" ~ "04",  # Arizona
    zip == "86315" ~ "04",  # Arizona
    zip == "94505" ~ "06",  # California
    TRUE ~ state  # keep existing value if zip doesn't match
  )
)





##### Step 3: baseline regressions #####



# setup regressions
base_controls <- 
  c(
  "pos_lin_msnbc",
  "i(chan_config)",
  "num_channels",
  "num_bc")

demog_controls <- 
  c(
    "pct_black", "pct_asian", "pct_other", "pct_hisp",
    "pct_male",
    "pct_age_10s", "pct_age_20s", "pct_age_30s", "pct_age_40s", "pct_age_50s", "pct_age_60s", "pct_age_70s", "pct_age_80s",
    "i(income_dec)",
    "pct_hsgrad", "pct_somecollege", "pct_bach", "pct_postgrad",
    "pct_suburban", "pct_urban",
    "totpop"
  )

for(dv in c(
  "amount_tea_party", "amount_non_tea_party", 
  "n_donors_tea_party", "n_donors_non_tea_party",
  "amount_dem", "n_dem"
)
) {
  
  print(paste0("-----", dv, "-----"))
  
  assign(paste0(dv, "_f_base"), fixest_form(y=dv, x=c("pos_lin_fnc", base_controls), fe="0"))
  assign(paste0(dv, "_f_controls"), fixest_form(y=dv, x=c("pos_lin_fnc", base_controls, demog_controls), fe="0"))
  assign(paste0(dv, "_f_controls_statefe"), fixest_form(y=dv, x=c("pos_lin_fnc", base_controls, demog_controls), fe="state"))
  
  assign(paste0(dv, "_m_base"), feols(get(paste0(dv, "_f_base")), data=chpos[position_year == 2009], cluster="stid"))
  assign(paste0(dv, "_m_controls"), feols(get(paste0(dv, "_f_controls")), data=chpos[position_year == 2009], cluster="stid"))
  assign(paste0(dv, "_m_controls_statefe"), feols(get(paste0(dv, "_f_controls_statefe")), data=chpos[position_year == 2009], cluster="stid"))
  
}



# Print Figure G.2.1 -----

resid_y <- resid(feols(fixest_form(y="amount_tea_party", 
                                   x=c(
                                     # "pos_lin_fnc",
                                     base_controls, demog_controls), fe="state"),
                 chpos[position_year == 2009], cluster="stid"))
resid_x <- resid(feols(fixest_form(y="pos_lin_fnc", 
                                   x=c(
                                     # "pos_lin_fnc",
                                     base_controls, demog_controls), fe="state"),
                       chpos[position_year == 2009], cluster="stid"))
X = "FNC Cable Channel Position"
Y = "Contribution Amount to Tea Party Candidates"
ggplot(data.frame(resid_y, resid_x), aes(x = resid_x, y = resid_y)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = "lm", se = FALSE, color = "blue") +
  labs(
    x = paste0("Residualized ", X),
    y = paste0("Residualized Log-Transformed\n", Y)
  ) +
  theme_light() +
  scale_y_continuous(trans = "log1p") +
  theme(    axis.text = element_text(size=12),
            axis.title = element_text(size=15),
            legend.title = element_text(size=15),
            legend.text = element_text(size=12)) %>% print()



# Estimate stacked regression 

temp <- chpos[position_year == 2009]
temp2 <- temp
temp$variable <- "tea_party"
temp2$variable <- "non_tea_party"
temp %<>% mutate(
  amount = amount_tea_party,
  n_donors = n_donors_tea_party
)
temp2 %<>% mutate(
  amount = amount_non_tea_party,
  n_donors = n_donors_non_tea_party
)
temp3 <- bind_rows(temp, temp2)
temp3 %<>% mutate(
  is_tea_party = as.integer(variable=="tea_party")
)
temp3 %<>% mutate(
  stacked_state = as.factor(paste(state, variable))
)

# setup stacked regressions
for(var in 
    c("pos_lin_fnc", 
      str_replace_all(str_replace_all(base_controls, "i\\(", ""), "\\)", ""), 
      str_replace_all(str_replace_all(demog_controls, "i\\(", ""), "\\)", "")
    )) {
  temp3[[paste0(var, "_tea")]] = temp3[[var]]*temp3[["is_tea_party"]]
}

stacked_base_controls <- 
  c(
    "pos_lin_msnbc", "pos_lin_msnbc_tea",
    "i(chan_config)", "i(chan_config_tea)",
    "num_channels", "num_channels_tea",
    "num_bc", "num_bc_tea")

stacked_demog_controls <- 
  c(
    "pct_black", "pct_black_tea", 
    "pct_asian", "pct_asian_tea", 
    "pct_other", "pct_other_tea", 
    "pct_hisp", "pct_hisp_tea",
    "pct_male", "pct_male_tea",
    "pct_age_10s", "pct_age_10s_tea", 
    "pct_age_20s", "pct_age_20s_tea", 
    "pct_age_30s", "pct_age_30s_tea", 
    "pct_age_40s", "pct_age_40s_tea", 
    "pct_age_50s", "pct_age_50s_tea", 
    "pct_age_60s", "pct_age_60s_tea", 
    "pct_age_70s", "pct_age_70s_tea", 
    "pct_age_80s", "pct_age_80s_tea", 
    "i(income_dec)", "i(income_dec_tea)",
    "pct_hsgrad", "pct_hsgrad_tea", 
    "pct_somecollege", "pct_somecollege_tea", 
    "pct_bach", "pct_bach_tea", 
    "pct_postgrad", "pct_postgrad_tea",
    "pct_suburban", "pct_suburban_tea", 
    "pct_urban", "pct_urban_tea",
    "totpop", "totpop_tea"
  )

for(dv in c("amount", "n_donors")) {
  
  print(paste0("-----", dv, "-----"))
  
  assign(paste0("stacked_", dv, "_f_base"), fixest_form(y=dv, x=c("is_tea_party", "pos_lin_fnc", "pos_lin_fnc_tea", stacked_base_controls), fe="0"))
  assign(paste0("stacked_", dv, "_f_controls"), fixest_form(y=dv, x=c("is_tea_party", "pos_lin_fnc", "pos_lin_fnc_tea", stacked_base_controls, stacked_demog_controls), fe="0"))
  assign(paste0("stacked_", dv, "_f_controls_statefe"), fixest_form(y=dv, x=c("is_tea_party", "pos_lin_fnc", "pos_lin_fnc_tea", stacked_base_controls, stacked_demog_controls), fe="stacked_state"))

  assign(paste0("stacked_", dv, "_m_base"), feols(get(paste0("stacked_", dv, "_f_base")), data=temp3[position_year == 2009], cluster="stid"))
  assign(paste0("stacked_", dv, "_m_controls"), feols(get(paste0("stacked_", dv, "_f_controls")), data=temp3[position_year == 2009], cluster="stid"))
  assign(paste0("stacked_", dv, "_m_controls_statefe"), feols(get(paste0("stacked_", dv, "_f_controls_statefe")), data=temp3[position_year == 2009], cluster="stid"))
  
}





##### Step 4: output baseline regressions #####



## setup for tables
setFixest_dict(c(zip="Zip",
                 state="State",
                 pos_lin_fnc = "FNC Channel Pos.",
                 pos_lin_msnbc = "MSNBC Channel Pos.",
                 amount_tea_party="Tea Party Candidates",
                 amount_non_tea_party="Other Rep. Candidates",
                 n_donors_tea_party="Tea Party Candidates",
                 n_donors_non_tea_party="Other Rep. Candidates",
                 amount_dem = "Dollar Amounts",
                 n_dem = "Number of Itemized Contributors",
                 stid = "Cable System"
),
 reset=TRUE
)


setFixest_etable(fitstat = ~ n + r2)

tablestyle  = style.tex(main="aer",
                        fixef.suffix = " FEs",
                        fixef.where ="var",
                        fixef.title = "",
                        stats.title = "\\midrule",
                        tablefoot=F,
                        yesNo="\\checkmark")



if(length(amount_tea_party_m_base$obs_selection$obsRemoved)==0) {
  mean_y1 <- mean(as.integer(chpos$amount_tea_party), na.rm = TRUE)
} else {
  mean_y1 <- mean(chpos$amount_tea_party[setdiff(1:nrow(chpos), -amount_tea_party_m_base$obs_selection$obsRemoved)], na.rm = TRUE)
}
mean_y2 <- mean(chpos$amount_tea_party[setdiff(1:nrow(chpos), -amount_tea_party_m_controls$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y3 <- mean(chpos$amount_tea_party[setdiff(1:nrow(chpos), -amount_tea_party_m_controls_statefe$obs_selection$obsRemoved)], na.rm = TRUE)
if(length(amount_non_tea_party_m_base$obs_selection$obsRemoved)==0) {
  mean_y4 <- mean(as.integer(chpos$amount_non_tea_party), na.rm = TRUE)
} else {
  mean_y4 <- mean(chpos$amount_non_tea_party[setdiff(1:nrow(chpos), -amount_non_tea_party_m_base$obs_selection$obsRemoved)], na.rm = TRUE)
}
mean_y5 <- mean(chpos$amount_non_tea_party[setdiff(1:nrow(chpos), -amount_non_tea_party_m_controls$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y6 <- mean(chpos$amount_non_tea_party[setdiff(1:nrow(chpos), -amount_non_tea_party_m_controls_statefe$obs_selection$obsRemoved)], na.rm = TRUE)

# Print Table 5 -----

print(etable(amount_tea_party_m_base, amount_tea_party_m_controls, amount_tea_party_m_controls_statefe,
       amount_non_tea_party_m_base, amount_non_tea_party_m_controls, amount_non_tea_party_m_controls_statefe,
       extralines = list("__Mean of DV" = c(sprintf("%.3f", mean_y1), sprintf("%.3f", mean_y2), sprintf("%.3f", mean_y3), sprintf("%.3f", mean_y4), sprintf("%.3f", mean_y5), sprintf("%.3f", mean_y6))),
       # file = "contributions_amount_tea_party_combined.tex",
       title = "FNC Effect on Zip Code-Level Total Itemized Contributions In Dollars",
       label = "table:contributions-amount-tea-party-combined",
       cluster=~stid,
       drop = "(Constant)",
       group=list("\\midrule Cable system controls"=c("chan", base_controls),
                  "Demographic controls"=c("income", demog_controls)),
       digits=3,
       digits.stats=2,
       style.tex=tablestyle,
       replace = T,
       depvar = T
))



if(length(n_donors_tea_party_m_base$obs_selection$obsRemoved)==0) {
  mean_y1 <- mean(as.integer(chpos$n_donors_tea_party), na.rm = TRUE)
} else {
  mean_y1 <- mean(chpos$n_donors_tea_party[setdiff(1:nrow(chpos), -n_donors_tea_party_m_base$obs_selection$obsRemoved)], na.rm = TRUE)
}
mean_y2 <- mean(chpos$n_donors_tea_party[setdiff(1:nrow(chpos), -n_donors_tea_party_m_controls$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y3 <- mean(chpos$n_donors_tea_party[setdiff(1:nrow(chpos), -n_donors_tea_party_m_controls_statefe$obs_selection$obsRemoved)], na.rm = TRUE)
if(length(n_donors_non_tea_party_m_base$obs_selection$obsRemoved)==0) {
  mean_y4 <- mean(as.integer(chpos$n_donors_non_tea_party), na.rm = TRUE)
} else {
  mean_y4 <- mean(chpos$n_donors_non_tea_party[setdiff(1:nrow(chpos), -n_donors_non_tea_party_m_base$obs_selection$obsRemoved)], na.rm = TRUE)
}
mean_y5 <- mean(chpos$n_donors_non_tea_party[setdiff(1:nrow(chpos), -n_donors_non_tea_party_m_controls$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y6 <- mean(chpos$n_donors_non_tea_party[setdiff(1:nrow(chpos), -n_donors_non_tea_party_m_controls_statefe$obs_selection$obsRemoved)], na.rm = TRUE)

# Print Table G.3.1 -----

print(etable(n_donors_tea_party_m_base, n_donors_tea_party_m_controls, n_donors_tea_party_m_controls_statefe,
       n_donors_non_tea_party_m_base, n_donors_non_tea_party_m_controls, n_donors_non_tea_party_m_controls_statefe,
       extralines = list("__Mean of DV" = c(sprintf("%.3f", mean_y1), sprintf("%.3f", mean_y2), sprintf("%.3f", mean_y3), sprintf("%.3f", mean_y4), sprintf("%.3f", mean_y5), sprintf("%.3f", mean_y6))),
       # file = "contributions_n_donors_tea_party_combined.tex",
       title = "FNC Effect on Zip Code-Level Total Number of Itemized Contributors",
       label = "table:contributions-n-donors-tea-party-combined",
       cluster=~stid,
       drop = "Constant",
       group=list("\\midrule Cable system controls"=c("chan", base_controls),
                  "Demographic controls"=c("income", demog_controls)),
       digits=3,
       digits.stats=2,
       style.tex=tablestyle,
       replace = T,
       depvar = T,
       placement = "H"
))



setFixest_dict(c(zip="Zip",
                 state="State",
                 pos_lin_fnc = "FNC Channel Pos.",
                 pos_lin_msnbc = "MSNBC Channel Pos.",
                 amount_tea_party="Tea Party Candidates",
                 amount_non_tea_party="Other Rep. Candidates",
                 n_donors_tea_party="Tea Party Candidates",
                 n_donors_non_tea_party="Other Rep. Candidates",
                 amount_dem = "Dollar Amounts",
                 n_dem = "Number of Itemized Contributors"
),
reset=TRUE
)



if(length(amount_dem_m_base$obs_selection$obsRemoved)==0) {
  mean_y1 <- mean(as.integer(chpos$amount_dem), na.rm = TRUE)
} else {
  mean_y1 <- mean(chpos$amount_dem[setdiff(1:nrow(chpos), -amount_dem_m_base$obs_selection$obsRemoved)], na.rm = TRUE)
}
mean_y2 <- mean(chpos$amount_dem[setdiff(1:nrow(chpos), -amount_dem_m_controls$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y3 <- mean(chpos$amount_dem[setdiff(1:nrow(chpos), -amount_dem_m_controls_statefe$obs_selection$obsRemoved)], na.rm = TRUE)

# Print Table G.5.1 -----

etable(amount_dem_m_base, amount_dem_m_controls, amount_dem_m_controls_statefe,
       extralines = list("__Mean of DV" = c(sprintf("%.3f", mean_y1), sprintf("%.3f", mean_y2), sprintf("%.3f", mean_y3))),
       # file = "contributions_amount_dem.tex",
       title = "FNC Effect on Zip Code-Level Total Itemized Contributions to Democratic Candidates",
       label = "table:contributions-amount-dem",
       cluster=~stid,
       drop = "Constant",
       group=list("\\midrule Cable system controls"=c("chan", base_controls),
                  "Demographic controls"=c("income", demog_controls)),
       digits=3,
       digits.stats=2,
       style.tex=tablestyle,
       replace = T,
       placement = "H"
) %>% print()



if(length(n_dem_m_base$obs_selection$obsRemoved)==0) {
  mean_y1 <- mean(as.integer(chpos$n_dem), na.rm = TRUE)
} else {
  mean_y1 <- mean(chpos$n_dem[setdiff(1:nrow(chpos), -n_dem_m_base$obs_selection$obsRemoved)], na.rm = TRUE)
}
mean_y2 <- mean(chpos$n_dem[setdiff(1:nrow(chpos), -n_dem_m_controls$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y3 <- mean(chpos$n_dem[setdiff(1:nrow(chpos), -n_dem_m_controls_statefe$obs_selection$obsRemoved)], na.rm = TRUE)

# Print Table G.5.2 -----

etable(n_dem_m_base, n_dem_m_controls, n_dem_m_controls_statefe,
       extralines = list("__Mean of DV" = c(sprintf("%.3f", mean_y1), sprintf("%.3f", mean_y2), sprintf("%.3f", mean_y3))),
       # file = "contributions_n_dem.tex",
       title = "FNC Effect on Zip Code-Level Total Number of Itemized Contributors to Democratic Candidates",
       label = "table:contributions-n-donors-dem",
       cluster=~stid,
       drop = "Constant",
       group=list("\\midrule Cable system controls"=c("chan", base_controls),
                  "Demographic controls"=c("income", demog_controls)),
       digits=3,
       digits.stats=2,
       style.tex=tablestyle,
       replace = T,
       placement = "H"
) %>% print()



setFixest_dict(c(zip="Zip",
                 state="State",
                 stacked_state="State-by-Tea Party",
                 is_tea_party = "\\multirow{2}{*}{Tea Party Candidates}",
                 pos_lin_fnc = "\\multirow{2}{*}{FNC Channel Pos.}",
                 pos_lin_fnc_tea = "\\multirow{2}{*}{FNC Channel Pos. $\\times$ Tea Party Candidates}",
                 pos_lin_msnbc = "\\multirow{2}{*}{MSNBC Channel Pos.}",
                 pos_lin_msnbc_tea = "\\multirow{2}{*}{MSNBC Channel Pos. $\\times$ Tea Party Candidates}",
                 amount_tea_party="Contributions to Tea Party Candidates ($)",
                 amount_non_tea_party="Contributions to non-Tea Party R Candidates ($)",
                 amount="Dollar Amounts",
                 n_donors_tea_party="No. Contributors to Tea Party Candidates",
                 n_donors_non_tea_party="No. Contributors to non-Tea Party R Candidates",
                 n_donors="Number of Itemized Contributors"
),
 reset=TRUE
)



tablestyle  = style.tex(main="aer",
                        fixef.suffix = " FEs",
                        fixef.where ="var",
                        fixef.title = "",
                        stats.title = "\\midrule",
                        tablefoot=F,
                        yesNo="\\checkmark")

set_col_widths = function(x, n, widths){
  # x: the character vector returned by etable
  # n: number of regs
  
  if(!missing(widths)){
    
    widths <- unlist(str_split(widths, ","))
    tex2replace <-
      str_replace(x, 
                 paste0("l", paste0(rep("c", (n-1)), collapse=""),""),
                 paste0("p{", widths[1], "cm}p{", 
                        paste0(widths[2:(n-1)], collapse="cm}p{"), 
                               "cm}p{",widths[n],
                               "cm}")
                 )
  }
  
    tex2replace
  
}



if(length(stacked_amount_m_base$obs_selection$obsRemoved)==0) {
  mean_y1 <- mean(as.integer(temp3$amount), na.rm = TRUE)
} else {
  mean_y1 <- mean(temp3$amount[setdiff(1:nrow(temp3), -stacked_amount_m_base$obs_selection$obsRemoved)], na.rm = TRUE)
}
mean_y2 <- mean(temp3$amount[setdiff(1:nrow(temp3), -stacked_amount_m_controls$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y3 <- mean(temp3$amount[setdiff(1:nrow(temp3), -stacked_amount_m_controls_statefe$obs_selection$obsRemoved)], na.rm = TRUE)
if(length(stacked_n_donors_m_base$obs_selection$obsRemoved)==0) {
  mean_y4 <- mean(as.integer(temp3$n_donors), na.rm = TRUE)
} else {
  mean_y4 <- mean(temp3$n_donors[setdiff(1:nrow(temp3), -stacked_n_donors_m_base$obs_selection$obsRemoved)], na.rm = TRUE)
}
mean_y5 <- mean(temp3$n_donors[setdiff(1:nrow(temp3), -stacked_n_donors_m_controls$obs_selection$obsRemoved)], na.rm = TRUE)
mean_y6 <- mean(temp3$n_donors[setdiff(1:nrow(temp3), -stacked_n_donors_m_controls_statefe$obs_selection$obsRemoved)], na.rm = TRUE)

# Print Table G.4.1 -----

etable(stacked_amount_m_base, stacked_amount_m_controls, stacked_amount_m_controls_statefe,
       stacked_n_donors_m_base, stacked_n_donors_m_controls, stacked_n_donors_m_controls_statefe,
       extralines = list("__Mean of DV" = c(sprintf("%.3f", mean_y1), sprintf("%.3f", mean_y2), sprintf("%.3f", mean_y3), sprintf("%.3f", mean_y4), sprintf("%.3f", mean_y5), sprintf("%.3f", mean_y6))),
       # file = "contributions_stacked_amount_combined.tex",
       title = "Stacked Regressions Comparing Zip Code-Level Total Itemized Contribution to Tea Party vs. Non-Tea Party Republican Candidates",
       label = "table:contributions-stacked-amount-combined",
       cluster=~stid,
       drop = "Constant",
       group=list("\\midrule Cable system-by-Tea Party controls"=c("chan", stacked_base_controls),
                 "Demographic-by-Tea Party controls"=c("income", stacked_demog_controls)),
       digits=3,
       digits.stats=2,
       style.tex=tablestyle,
       replace = T,
       placement="htbp"
) %>% print()




##### Step 5: sensitivity analysis #####



# setup regressions
base_controls <- 
  c(
    "pos_lin_msnbc",
    "i(chan_config)",
    "num_channels",
    "num_bc")


demog_controls <- 
  c(
    "pct_black", "pct_asian", "pct_other", "pct_hisp",
    "pct_male",
    "pct_age_10s", "pct_age_20s", "pct_age_30s", "pct_age_40s", "pct_age_50s", "pct_age_60s", "pct_age_70s", "pct_age_80s",
    "i(income_dec)",
    "pct_hsgrad", "pct_somecollege", "pct_bach", "pct_postgrad",
    "pct_suburban", "pct_urban",
    "totpop"
  )

for(dv in c(
  "amount_tea_party" 
)) {
  
  print(paste0("-----", dv, "-----"))
  
  # estimate reduced models on the same sample
  assign(paste0(dv, "_f_reduced"), fixest_form(y=dv, x=c("pos_lin_fnc", base_controls), fe="state"))
  
  full_model <- get(paste0(dv, "_m_controls_statefe"))
  rows.to.remove <- -full_model$obs_selection$obsRemoved

  assign(paste0(dv, "_m_reduced"), feols(get(paste0(dv, "_f_reduced")), 
                                         data=chpos[setdiff(1:nrow(chpos), rows.to.remove), ], cluster="stid"))
  
  # calculate oster ratio
  reduced_model <- get(paste0(dv, "_m_reduced"))
  coef_reduced <- reduced_model$coefficients["pos_lin_fnc"]
  r2_reduced <- as.numeric(fitstat(reduced_model, "wr2"))
  coef_full <- full_model$coefficients["pos_lin_fnc"]
  r2_full <- as.numeric(fitstat(full_model, "wr2"))
  oster_ratio <- coef_full*(r2_full-r2_reduced)/((coef_reduced-coef_full)*(min(1, 1.3*r2_full)-r2_full))
  assign(paste0(dv, "_oster_ratio"), oster_ratio)
  print(oster_ratio)
  # setwd(misc_dir)
  # write(round(get(paste0(dv, "_oster_ratio")), digits=1), paste0(dv, "_oster_ratio.tex"))
  
  # calculate oster bounds
  coef_bound <- coef_full - (coef_reduced-coef_full)*(min(1, 1.3*r2_full)-r2_full)/(r2_full-r2_reduced)
  print(coef_bound)
  stopifnot(coef_bound >= coef_full)
  if(round(coef_full, digits=1)==0) {
    oster_bound <- paste0("[", as.character(signif(coef_full, digits=3)), ", ", as.character(signif(coef_bound, digit=3)), "]")
  } else{
    oster_bound <- paste0("[", as.character(round(coef_full, digits=1)), ", ", as.character(round(coef_bound, digit=1)), "]")
  }
  print(oster_bound)
  # write(oster_bound, paste0(dv, "_oster_bound.tex"))
  
  ## setup for tables
  setFixest_dict(c(zip="Zip",
                   state="State",
                   pos_lin_fnc = "FNC Channel Pos.",
                   pos_lin_msnbc = "MSNBC Channel Pos.",
                   amount_tea_party="Tea Party Candidates",
                   amount_non_tea_party="Other Rep. Candidates",
                   n_donors_tea_party="Tea Party Candidates",
                   n_donors_non_tea_party="Other Rep. Candidates",
                   amount_dem = "Dollar Amounts",
                   n_dem = "Number of Itemized Contributors",
                   stid = "Cable System"
  ),
  reset=TRUE
  )
  
  
  setFixest_etable(fitstat = ~ n + r2)
  
  tablestyle  = style.tex(main="aer",
                          fixef.suffix = " FEs",
                          fixef.where ="var",
                          fixef.title = "",
                          stats.title = "\\midrule",
                          tablefoot=F,
                          yesNo="\\checkmark")

  mean_y1 <- mean(chpos[[dv]][setdiff(1:nrow(chpos), -full_model$obs_selection$obsRemoved)], na.rm = TRUE)
  mean_y2 <- mean(chpos[[dv]][setdiff(1:nrow(chpos), -full_model$obs_selection$obsRemoved)], na.rm = TRUE)
  
  # Print Table G.1.1 -----
  
  print(etable(reduced_model, full_model,
         extralines = list("__Mean of DV" = c(sprintf("%.3f", mean_y1), sprintf("%.3f", mean_y2))),
         # file = paste0("contributions_", ifelse(dv=="amount_tea_party", "amount_tea_party", "n_donors_tea_party"), "_oster_ratio.tex"),
         title = paste0("FNC Effect on Zip Code-Level Total ", ifelse(dv=="amount_tea_party", "Itemized Contributions In Dollars", " Number of Itemized Contributors"), " (for Oster Ratio Calculation)"),
         label = paste0("table:contributions-", ifelse(dv=="amount_tea_party", "amount-tea-party", "n-donors-tea-party"), "-oster-ratio"),
         cluster=~stid,
         drop = "(Constant)",
         group=list("\\midrule Cable system controls"=c("chan", base_controls),
                    "Demographic controls"=c("income", demog_controls)),
         digits=3,
         digits.stats=2,
         style.tex=tablestyle,
         replace = T,
         depvar = T,
         placement = "H"
  ))

}


